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ABSTRACT 


The  use  of  the  Fisher  distribution  to  estimate 

the  mean  direction  of  magnetization  of  a  rock  at  a 

sampling  site  is  now  standard.  Sampling  sites  are 
4  5 

chosen  to  cover  10  to  10  years  to  average  out  the 
effect  of  secular  variation.  The  controversy  about 
how  to  combine  these  site  means  has  never  been  satis¬ 
factorily  resolved.  By  using  statistical  models  for 
secular  variation,  this  paper  suggests  how  methods 
should  be  derived.  A  number  of  interesting  statis¬ 
tical  distributions  and  estimation  problems  are 
shown. 
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1.  INTRODUCTION 

When  a  lava  containing  iron  cools  below  its  Curie  point,  it  becomes 
magnetized  in  the  direction  of  the  earth's  field  at  that  place  and  time. 
Sediments  containing  magnetic  particles  also  acquire  (a  much  weaker)  local 
magnetization  as  they  are  formed.  The  study  of  palaeomagnetism  led  to  the 
current  revolution  in  geology  -  plate  tectonics.  For  assuming  that  the  earth's 
magnetic  field  has  always  been  much  the  same  as  it  is  now  -  except  for  rever¬ 
sals  of  polarity  -  a  magnetic  measurement  gives  a  good  estimate  of  the  lati¬ 
tude  of  the  site  when  the  rock  was  magnetized.  These  latitudes  make  no  sense 
unless  one  is  willing  to  admit  that  continental  drift  has  taken  place  -  in 
fact,  they  allow  continental  reconstructions  to  be  made.  A  full  account  may 
be  found  in  the  book  by  McEl hinny  (1973). 

The  strength  of  the  magnetization  of  a  specimen  may  have  altered  over 
time  and  so  it  is  not  used  in  this  work  -  only  the  direction  of  natural  reman¬ 
ent  magnetism,  i.e.,  the  original  magnetization.  Later  changes  can  be 
"cleaned"  out  by  methods  not  relevant  here.  (If  the  rocks  are  unstably  mag¬ 
netized,  they  cannot  be  used.)  If  N  specimens  are  taken  from  one  site 

with  directions,  the  N  unit  vectors,  r^ .  r^,  after  cleaning,  may  be 

assumed  to  be  a  sample  from  Fisher's  distribution 

4*  s1nh~  e*p  K  r'p  <K  2  °>  (I) 

where  u  Is  the  (unit)  mean  direction.  Fisher  (1953)  showed  that  the  maximum 
likelihood  estimator  of  p  is  the  direction  of  Re  £r*,  the  vector  resultant 
of  the  sample  and,  further,  that  if  |R|  is  the  length  of  R,  ic  is  the  solu- 


A 

An  approximation  to  <  is  <  given  by 


< 


N  -  1 

N  -  P*T  * 


(2  ) 


k  measures  the  precision  of  the  Fisher  distribution.  Further  statistical 
methods  appropriate  for  palaeomagneticists  were  given  by  Watson  (1956a,  b) 
and  Watson  and  Irving  (1957). 

The  latter  paper  gives  an  approximate  within  and  between  sites  method 

of  analyzing  data  from  several  sites.  It  is  assumed  that  the  within  sites 

distribution  is  Fisher  with  <  -  k  about  the  site  mean  and  that  the  site 

w 

means  are  Independently  Fisher  with  k  *  <g  about  the  true  palaeomagnetic 

direction.  Unfortunately,  this  compound  of  Fisher  distributions  is  only 

approximately  Fisher.  The  estimate  of  Kg  was  meant  to  measure  the  secular 

variation  since  it  is  supposed  that  the  sites  sample  the  full  cycle  of  secu- 

4  5 

lar  variations.  This  is  some  10  to  10  years  when  the  pole  moves  around 
its  mean  position,  a  short  period  on  continental  drift  time  scales. 

In  designing  this  analysis,  no  real  thought  was  given  to  modelling  secu¬ 
lar  variation.  Furthermore,  the  immense  amount  of  experimental  work  done 
since  then  has 'shown  that  site  means  often  do  not  vary  symmetrically  around 
some  mean  direction  as  implied  by  the  Fisher  distribution  (1),  as  we  supposed. 

When  the  earth's  field  is  averaged  over  the  time  period  of  secular  var¬ 
iation,  It  Is  found  to  be  approximated,  with  surprising  accuracy,  by  an  earth- 
centered  magnetic  dipole  currently  Inclined  at  about  11  degrees  to  the  rota¬ 
tional  axis.  From  elementary  physics,  the  magnetic  force  H  at  a  point  r 
from  a  dipole  M  (H,  M,  r  are  all  vectors)  is  given  by 


H 


(3) 


where  the  prime  denotes  a  transpose,  I  is  the  3x3  identity  matrix,  and 


| r |  is  the  length  of  r.  Of  course,  the  center  of  the  earth  is  too  hot  to 
support  a  dipole!  The  field  is  thought  to  be  due  to  a  self-exciting  dynamo 
made  of  asymmetric  flows  of  conducting  rock  (mainly  iron  and  nickel)  through 
stray  permanent  fields  of  the  rotating  earth.  Minor  changes  in  the  flows  and 
the  eastward  movement  of  the  crust  relative  to  the  interior  (core)  of  the 
earth  cause  the  secular  variation  and  even  the  reversals. 

Irving  and  Ward  (1964)  supposed  that  the  secular  variation  was  due  to 
a  smaller  geocentric  dipole  m  whose  orientation  varied  uniformly  at  random 
over  the  period  but  whose  magnitude  |m|  was  constant.  Then  the  field  at 
site  r  at  different  times  in  the  period  is  independent  and  has  the  repre¬ 
sentation 

H  *  f3  "  *](*  +  "O/M3  (4) 

l  kr  J 

Here  and  below  we  will  regard  M  as  fixed  and  define  y  *  M/|m|.  Creer  et 
al.  (1959)  assume  that  the  main  dipole  M  has  fixed  strength  |M|  but 
“wobbles”  -  that  the  direction  of  the  main  dipole  is  independent  from  time  to 
time  and  has  a  Fisher  distribution.  Creer  et  al.'s  model  says 


where  M*  *  |M|d  and  d  has  a  Fisher  distribution  about  y,  and  some  large 
k.  Cox  (1970)  supposes  that  the  main  dipole  oscillates  In  strength  and  wobbles 
In  direction  and  that  there  Is,  In  addition  and  Independently,  a  randomly 
oriented  smaller  dipole  m,  as  in  the  Irving  and  Watson  model.  In  his  case 


-  i]m*  / 1 r | 3  (5) 

M 


3^-1 


M 


(M+  +  m)/ |r | 3 


(6) 


where  the  average  value  of  M+  Is  to  be  M  about  which  M+  has  a  rotationally 
symmetric  distribution. 
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In  the  last  five  years  some  more  complex  models  have  been  suggested 
which  generate  secular  variations  with  dipoles  on  the  surface  of  the  core  - 
see,  e.g.,  Harrison  &  Watkins  (1979).  Here  we  will  be  content  to  explore 

models  expressed  In  (4),  (5),  and  (6),  or  closely  related  to  them.  We  will 

assume  that  large  samples  have  been  taken  at  each  site  so  that  site  mean  direc 
tlons  are  known  exactly.  Let  these  directions  be  Lj,  Lg,  ....  LN-  The 
problem  is  to  estimate  y,  the  direction  or  mean  direction  of  M,  M*.  M+. 

We  will  choose,  as  a  unit  of  length,  the  earth's  radius  (assuming  it  to 
be  a  sphere)  and  set  u  «  r/|r|.  Then  the  common  factor  in  (4),  (5),  and  (6) 
is  3uu'  -  I  where  u  Is  a  fixed  but  arbitrary  unit  vector.  We  will  write 

U  »  3mT  -  I  .  (7) 

Note  that  U  *  U’,  *  3uu'  +  I  and  that 

U  *  2uu*  +  UgUg  +  u3Uj 

where  u,  u^,  u^  are  orthonormal.  Hence  the  eigenvalues  and  eigenvectors  of 
U  are  trivially  known.  Each  of  the  models  has  the  form 

H  «  UX 

where  X  is  a  random  vector,  U  is  fixed  and  only  the  direction  of  H  is 
observed.  In  the  next  section  we  will  explore  the  distributional  and  estima¬ 
tion  problems  this  raises.  They  have  some  Intrinsic  statistical  interest 
and  may  have  other  applications.  In  the  third  section,  we  return  explicitly 
to  our  main  problem. 

2.  STATISTICAL  PRELIMINARIES 

(a)  Fisher  observed  that  his  distribution  may  be  derived  from  the 

Gaussian  as  follows.  If  X  has  a  trlvarlate  Gaussian  distribution  with 

o 

mean  vector  y  and  covariance  matrix  o  I,  set  R  ■  |X|,  L  ■  X/R, 


X  *  y/|y|.  Then  the  joint  density  of  R  and  L  Is 


R  _  1 


(R2  + 


’)  exp  ^4^-  L'X 

G 


(8) 


so  that  the  density  of  L  on  the  unit  sphere,  conditional  on  a  fixed  R,  Is 
proportional  to 

exp-*4^l'X  .  (9) 

G 

From  (1)  we  see  that  L  then  has  a  Fisher  distribution  about  the  mean  direction 

2 

X  with  a  k  of  R|y|/o  .  If  (9)  were  appropriate,  we  know  how  to  estimate 
2 

X  and  R|y|/o  -  use  Fisher's  estimations  for  (1). 

If,  however,  (8)  were  appropriate  but  we  were  only  given  the  directions 
Lj,...,  Ln  of  Xj,...,  Xjy,  we  would  have  to  use  the  density  f  of  L, 

r 

2 

- T7T"T  exP - V  O'  +  *  2R  |  y  |  L'X)dR  .  (10) 

J  (2*)  '  o  2<T 
0 

This  Is  simpler  If  we  set  S  *  R/o,  v  ■  |y|/o,  since  then 


f(L.X,v) 


>00 

2 

— ^75-  exp{-  1(S2  +  V2  -  2SvL'X)JdS 

0 


with  the  vector  derivative 


if 

3X 


S3v 


(2ir) 


37? 


(u) 


(12) 


«  vJ(L,X,v)L,  say  .  (12') 

Since  X'X  •  1,  the  maximum  likelihood  (m.1.)  estimate  of  X  requires  us  to 
solve 


fr 
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Jr  (  Z  log  f(L .)  +  ex'x)  -  0 
dA  j*l  J 


l.e. , 


N  J(L,,A,v) 


■  (L<j  »v)  j  v 


l,  ♦22i-o 


If  v  is  known,  (13)  may  be  solved  numerically  by  an  iteration.  An  initial 
(1)  N 

estimate  of  X,  X  would  be  the  direction  of  I  L.  With  computer  programs 

1  3 

to  compute  J(Lj,A'^)  and  f(Lj,A^),  (13)  would  be  used  to  find  A^',  etc. 
The  two  integrals  have  well-behaved  integrands  so  that  numerical  integration 
will  not  be  difficult. 

To  estimate  the  non-negative  parameter  v,  we  observe  that 

xto 

2 

&  *  ^377  (-  V  ♦  Sl'X)  e*p{-  fa*  *  V2  -  2SvL'X)|dS 

J 


-vf  +  L'AO  . 


Hence  the  m.l.  equation  for  v  is 

N  vf(L. ,A,v)  +  La  J(L,,A,v) 

I - “ - m — r*- -T - J - *  0 

j»l  PTlJXS] 


.  N  L.'A  J(L.,X,v) 

±  J  ~  J.  , - \>  fit 

N  j.j  f(Lj.X,v)  v  *  U* 

a  form  Ideal  for  Iteration.  It  is  convenient  that  no  more  integrals  need  be 
evaluated.  To  start  off  the  v  iteration  we  again  use  the  analogy  between 
(8)  and  (9)  and  suggest 

v(,)  ■  fHtjr  •  <>' 
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Thus  to  solve  jointly  (13)  and  (15)  for  X  and  v,  one  seesaws  between 

them  using  the  suggested  X'^  and  v'^.  Hence  given  only  the  directions 

o 

Lj,...,  L|y  of  a  sample  of  N  from  G(y,o  I),  we  may  estimate  the  direction 

of  y  and  v  =  |y|/a.  To  find  the  covariance  matrix,  the  second  derivatives 
of  the  log-likelihood  will  be  evaluated  in  the  usual  way.  We  will  leave  this 
calculation  for  an  applied  paper  to  appear  elsewhere.  This  paper  will  also 
show  the  shape  of  distribution  defined  by  (10)  or  (11)  and  compare  it  with  others. 
This  distribution  has  of  course  appeared  before  (see  e.g.  Kendall  (1974)  under 
the  names  "off-set"  or  "displaced"  normal  —  I  prefer  "angular"  normal. 

(b)  If  we  observed  the  directions  of  n  copies  of  Y  *  TX,  where  X  is 
Gaussian  mean  y  with  known  covariance  matrix  o' I,  and  T  is  a  known  non- 

p 

singular  3*3  matrix,  we  see  that  we  may  set  o  =  1  without  loss  of  gen¬ 
erality  where  we  wish  to  estimate  the  direction  of  y.  The  density  of  Y  is 


(2it)3/2|  |T| | 


exp  -  y(Y'(TT')'1¥  -  2p'T'lV  +  u'u)  .  (17) 


Setting  Y  s  RL  where  R  >  0  and  |L|  *  1,  the  joint  density  of  R  and  L  is 

n2 


(2x)3/2||T|| 

Thus  the  density  of  L  is 


exp  -  i(R2L'(TT')"1L  -  2R.T_1L  +  u'u) 


-1, 


h(L,y) 


- 1™ - exp  -  ■5-(R2L'(TTT1L  -  2RyT~*L  +  y'y)dR 

(2tt)3/Z|  |T|  |  2 


so 


8h 

3V 


(2ir) 


3/2njjJ  *XP  ‘  ‘  2Ry'T"1L  +  y'y) 


( RT~1L  -  y)dR  . 
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Define 


k(L,u) 


3 

- JL - exp  -  •kR2L'(TTT1L  -  +  u'ujdR  . 

(2ir)3/Z||T||  21 

0 


We  see  that  the  m.l.  equation  for  y  is 


N 

z 


<(Lj,u)  T-l, 

Cm  T  L 


»  ^  HTtpr 


(18) 


A  natural  initial  value  is  y^  *  N“*T”*ZLj  for  the  iteration  to  find  y. 
Then  we  will  take  the  direction  of  y.  To  get  the  accuracy  of  y,  we  may  use 
second  derivatives  in  the  usual  way. 


(c)  To  go  one  step  further,  if  Y  *  UZ  where  Z  is  Gaussian  with  mean 
vector  y  and  a  known  covariance  matrix  Z  and  we  wish  to  estimate  the 
direction  of  y  given  the  direction  of  N  Y's,  we  may  write 


Y  *  UZ  *  UZ+1/2Z"1/2 
*  TX 


Z 


where 


T  »  UZ1/2  ,  X  =  I'1/2Z 


and  X  is  Gaussian  (Z~1/,2y,I).  Thus  we  may  use  the  method  of  subsection  (b) 
to  estimate  Z~^2y.  Given  an  estimate  of  the  latter,  we  have  an  estimate  of 
y,  since  Z^2  is  known,  and  hence  of  its  direction. 

Subsection  (c)  leads  to  an  algorithm  to  solve  the  problem  of  Section  1 
for  a  model  of  secular  variation  where  the  geocentric  dipole  has  a  multivariate 
Gaussian  distribution.  This  Is  not  a  model  Included  In  (4)  and  (5);  it  Is  a 
limiting  case  in  (6). 
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(d)  To  deal  with  model  (6)  we  are  led  to  study  the  distribution  of 
Y  *  X  +  v  where  X  is  Gaussian  -  mean  u»  covariance  l  -  and  v  is  uni¬ 
formly  distributed  over  a  sphere  of  radius  6.  This  leads  to  a  plethora  of 
apparently  new  multivariate  densities. 

Because  of  their  intrinsic  interest,  we  consider  some  veAy  spec ajoJL  oaaea 
not  directly  relevant  to  our  main  problem.  In  one  dimension,  the  density  of 
y  clearly  is 

|-(exp  -  j(y  -  u  -  <5)2  +  exp  -  |<y  -  y  +  6)2) 

VXTT 


=  expj-  |c(y  -  u)2  +  62]|cosh  5(y  -  y)  . 


In  two  dimensions  the  density  of  y  *  (yj,y2K  is 


f2lT 

1  2  1  2 
exp  -  ^(y^  -  pj  -  8  cos  9)  -  ^(y2  -  u2  -  5  sin  9)  de 


^  exp  -  -^(y,  -  u,)2  +  (y9  -  u9)2  +  S2} 


i  t  v^2 


*  57 


■2it 

exp  {^(yj  -  ux)  cos  9  +  6(y2  -  u2)  sin  ejd9 


*  exp{“  ?  +  ^2  ’  v2 ^  + 

where  Iq(z)  is  the  imaginary  Bessel  function  of  zero  order. 

In  three  dimensions,  a  similar  averaging  of  the  Gaussian  with  mean  y 
and  covariance  matrix  I  leads  to 
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,xp  {.  fa  .  uny  .  „  .  *  Ml  . 

In  each  case,  as  6  -*•  0  the  modifying  factor  tends  to  unity.  The 
averaging  has  most  effect  when  |y  -  y|  for  the  basic  density  is  then  fallen 
off  very  fast.  Let  us  now  return  to  our  real  problem. 

Model  (6)  requires  a  general  form  of  the  interesting  new  distributions 
just  given.  To  put  it  in  neutral  notation,  let  Y  =  X  +  d  where  d  is  uni¬ 
formly  distributed  over  a  sphere  of  radius  6,  X  is  Gaussian  mean  y  and 
covariance  £,  and  the  density  of  Y  at  y  is  given  by 


Ave 

|d|-« 


exp  -  |<y  -  d  -  uKl_1(y  -  d  -  y) 


*  — Kj*  —  exp  -  i(y  -  y)z_1(y  -  y)  (19) 

(2tt)3/z  | Z |  2 


*|d|*6  exP(d^_1(y  -  P)  -  \  d'E-1d}  . 

If  E  *  H'DH  where  H  is  orthogonal  and  D  diagonal,  z  =  6D“*H(y  -  y),  and 
-1  2  -1 

E  s  5  D  ,  the  averaging  term  is 

iJflj  exp  {e'z  •  j  e'E^e), 

which  cannot  be  further  evaluated.  In  fact,  with  z  =  0,  it  is  the  normal¬ 
izing  constant  of  the  Bingham  distribution  (see,  e.g.,  the  book  by  Mardia,  1972) 
and,  when  z  t  0,  of  a  generalization  recently  studied  by  Beran  (1979).  To 
get  to  the  working  form  of  the  Cox  model,  we  must  find  the  density  of  the 
direction  of  the  vector  y  from  (19).  Thus  the  Cox  model  (6)  Is  hard  to 
deal  with  mathematically,  so  estimation  schemes  based  on  it  do  not  seem  very 


practical.  However,  the  Gaussian  model  without  the  random  smaller  dipole 
seems  an  adequate  approximation  since,  typically,  |m|/|M(  is  about  one-tenth. 

(e)  The  model  (4)  raises  interesting  problems.  Here  M  is  fixed  and 
m  is  uniformly  distributed  over  a  sphere  of  radius  |m|,  also  fixed.  Thus 
Em  =  0,  Emm"  =  |m|l/3. 

First,  we  consider  several  related  academic  problems  similar  to  familiar 
textbook  estimation  examples.  Suppose  that  N  >  2  points  y  are  uniformly 
distributed  along  a  line  segment  in  space,  i.e.,  y  =  M  +  vd  where  v  is 
uniformly  distributed  on  (0,|m|)  with  M,  |m|  and  d  unknown  and  to  be 
estimated.  If,  when  the  points  are  arranged  on  the  segment. they  fall  in  the 
order  yj,yn,...,  yN,  then  +  d  is  known  exactly,  and  we  have 
ly2  “  yll’ly3  "  yz\*  •••»  '  yN  I  as  the  intenor  9aPs  when  N  points 

are  randomly  distributed  on  (0,|m|).  The  sufficient  statistic  is  (y^y^)- 
Simple  estimators  are  (with  m  *  d|m|)  y^  =  M,  y^  »  M  +  m  or  y^  *  M  -  m, 
y^  *  M.  These  will  be  biased,  as  are  the  analogues  on  the  real  line.  To  give 
another  interesting  example,  suppose  we  could  observe'  N  copies  of  y  =  M  +  m 
where  M  and  m  are  as  in  model  (4).  For  N  z  3,  |m|  and  M  can  be  found 
exac tty  from  the  perpendicular  bisectors  of  chords  which  must  intersect  at 
M.  As  is  so  often  the  case,  things  are  simpler  on  the  circle  and  sphere  than 
on  the  line. 

In  practice,  we  can  only  observe  the  directions  L  of  U(M  +  m).  Assuming 
that  M  *  |m|p/c,  m*  |m|d,  with  c  ■  |m|/|M|  known  and  cpproximately  0.1  and 
d  uniformly  distributed  over  the  surface  of  the  unit  sphere.  It  suffices  to 
consider  the  case  of  observing  L  where 

RL  -  U(u  +  cd) 

with  known  c.  M.l.  here  is  very  awkward.  Since 
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y  +  cd.  -  RiU"1Li  (i  =  1,...,  N), 

one  might  use  (since  Ed  *  0)  the  estimator 

~  N  -i 

y  «  I  U  ll.  .  (20) 

1  1 

A  case  of  confidence  can  be  obtained  from  the  approximate  multivariate  dis- 
N  _i 

tribution  of  I  U  L. .  To  this  end  we  note  that 
1  1 

R2  *  y'U2y  +  2cy'U2d  +  c2d'U2d 


SO 


ER2  *  3(y'u)2  +  1  +  2c2 


where  2cr  may  be  ignored.  If  u  is  an  approximate  pole  position,  define 

R2  *  3({j'u)2  +  1, 

and  write 

y  +  cd  *  R  U_1L  . 

Then 

y  »  R  E  U_1L 
R  U_1L  -  y  *  cd 


SO 

E(R  U_1L  -  y)(R  U_1L  -  v)'  •  (c2/3)I  (21) 

which  leads  easily  to  the  procedure. 


(f)  The  model  (5)  of  Creer  et  al.  has  M  ■  |M|d  where  |M|  is  fixed 
and  known*  and  d  has  a  Fisher  distribution  about  y  with  accuracy  <c.  If 
the  direction  of  H  *  UM  is  known,  l.e.,  the  direction  L  of  Ud  is  known, 
how  to  estimate  y?  Since  ■  Ud^,  d^  »  R^lf1!^,  we  may  compute  each 

reduce  to  unit  length  (so  not  knowing  R  does  not  matter)  to  obtain  the 
dj.  Since  these  have  a  Fisher  distribution,  Fisher's  estimator  and  method  of 
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getting  a  confidence  case  for  p  may  be  used.  This  is,  in  fact,  a  standard 
method  of  estimating  the  positon  of  the  pole.  It  is  important  to  observe  that 
the  Fisherian  method  is  not  applied  here  to  the  site  means  but  their  transforms. 
Note  that  It  Is  related  to  but  different  from  estimators  (20),  (18),  (15),  (13), 
which  are  very  similar. 

3.  ESTIMATION  OF  POLE  POSITIONS 

We  have  seen  that  a  rational  method  of  estimating  pole  positions  should 
follow,  by  the  application  of  maximum  likelihood,  from  a  statistical  model 
for  secular  variation.  In  cases  where  the  model  makes  strict  m.l.  estimation 
too  difficult,  one  must  try  to  improve  and  check  the  approximations  and  simpli¬ 
fications  made.  Since,  however,  there  is  unlikely  to  be  any  agreement  on  the 
model,  the  chosen  method  should  not  be  sensitive  to  model  changes  within  the 
range  of  disagreement,  i.e.,  the  method  should  be  robust. 

The  multivariate  normal  method  of  Section  2(c)  allows  a  wide  variety  of 
models  to  be  tried.  These  will  be  tried  out  in  an  applied  paper  to  see  what 
model  changes  are  most  crucial.  The  Fisher  estimation  method  has  been  shown 
to  be  robust  against  deviations  from  the  Fisher  distribution  but  the  cone  of 
confidence  is  not  so  reliable  -  Watson  (1967). 

A  counsel  of  greater  perfection  is  to  make  a  more  detailed  statistical 
study  of  secular  variation  in  the  hope  of  finding  a  more  convincing  model. 
Certainly  the  newer,  more  complex  models  mentioned  earlier  should  be  explored. 
This  will  be  attempted  elsewhere.  Further,  we  have  here  assumed  that  the 
site  means  are  known,  whereas  they  will  only  be  estimated,  often  with  different 
accuracies.  This  must  be  taken  into  account  in  the  final  practical  method. 


1 
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